
setThreadOptions(numThreads=20,stackSize=1e07)

maindir<-getwd()

for(modeltype in c("simplified"#,"rich"
                   )){

  setwd(maindir)
  
source("./auxfiles/loaddata.R")

setwd(maindir)

varcov<-as.data.table(read.table("./output/varcovmat_simplified.txt"))

varcov<-varcov[!colnames(varcov)%in%c("hwreg_mean_married","cons_mean_married","cons_mean_age","cons_mean_age2"),]
varcov<-varcov[,!colnames(varcov)%in%c("hwreg_mean_married","cons_mean_married","cons_mean_age","cons_mean_age2"),with=FALSE]
testlev<-1


names(upper)<-names(newguess)
names(lower)<-names(newguess)

newguess[which(newguess>upper)]<-upper[which(newguess>upper)]
newguess[which(newguess<lower)]<-lower[which(newguess<lower)]

#Shutting down PT work moments:
newguess[grepl("part",names(newguess))]<-0
upper[grepl("part",names(newguess))]<-0
lower[grepl("part",names(newguess))]<-0


paramnums<-which(upper!=lower)

pos<-1
Jmat<-NULL
Cfactmat_ss <- NULL
Cfactmat_fs <- NULL
Cfactmat_stringency <- NULL
Cfactmat_marstringency <- NULL
Cfactmat_singstringency <- NULL
Cfactmat_earlydi <- NULL
Cfactmat_meanstest <- NULL
Cfactmat_meanstest20k <- NULL
Cfactmat_meanstest40k <- NULL
Cfactmat_meanstest90k <- NULL
Cfactmat_ss_meanstest <- NULL
Cfactmat_fs_meanstest <- NULL
Cfactmat_stringency_meanstest <- NULL
Cfactmat_marstringency_meanstest <- NULL
Cfactmat_singstringency_meanstest <- NULL
Cfactmat_earlydi_meanstest <- NULL
Cfactmat_fastdi <- NULL
Cfactmat_fastdi_meanstest <- NULL

Cfactlevels_ss <- NULL
Cfactlevels_fs <- NULL
Cfactlevels_stringency <- NULL
Cfactlevels_marstringency <- NULL
Cfactlevels_singstringency <- NULL
Cfactlevels_earlydi <- NULL
Cfactlevels_meanstest <- NULL
Cfactlevels_meanstest20k <- NULL
Cfactlevels_meanstest40k <- NULL
Cfactlevels_meanstest90k <- NULL
Cfactlevels_ss_meanstest <- NULL
Cfactlevels_fs_meanstest <- NULL
Cfactlevels_stringency_meanstest <- NULL
Cfactlevels_marstringency_meanstest <- NULL
Cfactlevels_singstringency_meanstest <- NULL
Cfactlevels_earlydi_meanstest <- NULL
Cfactlevels_fastdi <- NULL
Cfactlevels_fastdi_meanstest <- NULL

#Size of policy adjustment for counterfactuals
marginval <- 1.1

for(n in paramnums){
  if(names(paramnums[paramnums==n]) %in% c("age","spage")) boundnum <- 0.0025
  else if(names(paramnums[paramnums==n]) %in% c("age2","spage2")) boundnum <- 0.01
  else if(names(paramnums[paramnums==n]) %in% c("wagevar","spwagevar")) boundnum <- 0.00025
  else if(names(paramnums[paramnums==n]) %in% c("wagecorr")) boundnum <- 0.015
  else if(names(paramnums[paramnums==n]) %in% c("noisevar","spnoisevar")) boundnum <- 0.001
  else boundnum <- 0.025
  param1<-min(newguess[n]+boundnum,upper[n])
  param0<-max(newguess[n]-boundnum,lower[n])
  newguess1<- newguess
  newguess1[n]<-param1
  newguess0<- newguess
  newguess0[n]<-param0
  
  ###### Upper bound #########
  temp<-wagesetup_fit(newguess1,0,0,wage_gridsize=10,spwage_gridsize=5#,filesuffix=typename
  )
  #Wage<-temp$Wage
  #spWage<-temp$spWage
  wagepars<-temp$wagepars
  spwagepars<-temp$spwagepars
  rm(temp)
  
  temp<-paramsetup(newguess1,wagepars,spwagepars,
             ptwork_halfcost=ptwork_halfcost,
             single_samecost=single_samecost,
             ptwork_halfdis=ptwork_halfdis,
             single_samedis=single_samedis,
             spouse_samejobrisk=spouse_samejobrisk,
             single_sameprefs = single_sameprefs,
             health_simpleprefs = health_simpleprefs,
             single_sameallowanceprob=single_sameallowanceprob,
             work_noallowanceprob=work_noallowanceprob)
  params<-temp$params
  param_input<-temp$param_input
  rm(temp)
  shocks.wagenoise <- rnorm(50*sum(numsims),0,wagepars$noisevar)
  shocks.spwagenoise <- rnorm(50*sum(numsims),0,spwagepars$noisevar)
  shocks.wage <- exp(shocks.wageraw*wagepars$wagevar)
  shocks.spwage <- exp((shocks.spwageraw*(1-wagepars$wagecorr^2)^0.5 + shocks.wageraw*wagepars$wagecorr)*spwagepars$wagevar)
  
  
    simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  moments1<-const_moments(simdata, shocks.wagenoise=shocks.wagenoise,
                          shocks.spwagenoise=shocks.spwagenoise,
                          nocomp=nocomp,noeventDI=noeventDI,onlyFT=onlyFT,noflow=noflow,nostockDI = nostockDI, nocontflow = nocontflow,
                          consmoments_simple=moments_simple,
                          workmoments_simple=moments_simple,
                          DImoments_simple=moments_simple)
  
  ###### Lower bound #########
  temp<-wagesetup_fit(newguess0,0,0,wage_gridsize=10,spwage_gridsize=5#,filesuffix=typename
  )
  #Wage<-temp$Wage
  #spWage<-temp$spWage
  wagepars<-temp$wagepars
  spwagepars<-temp$spwagepars
  rm(temp)
  
  temp<-paramsetup(newguess0,wagepars,spwagepars,
                   ptwork_halfcost=ptwork_halfcost,
                   single_samecost=single_samecost,
                   ptwork_halfdis=ptwork_halfdis,
                   single_samedis=single_samedis,
                   spouse_samejobrisk=spouse_samejobrisk,
                   single_sameprefs = single_sameprefs,
                   health_simpleprefs = health_simpleprefs,
                   single_sameallowanceprob=single_sameallowanceprob,
                   work_noallowanceprob=work_noallowanceprob)
  params<-temp$params
  param_input<-temp$param_input
  rm(temp)
  shocks.wagenoise <- rnorm(50*sum(numsims),0,wagepars$noisevar)
  shocks.spwagenoise <- rnorm(50*sum(numsims),0,spwagepars$noisevar)
  shocks.wage <- exp(shocks.wageraw*wagepars$wagevar)
  shocks.spwage <- exp((shocks.spwageraw*(1-wagepars$wagecorr^2)^0.5 + shocks.wageraw*wagepars$wagecorr)*spwagepars$wagevar)
  
  simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  moments0<-const_moments(simdata, shocks.wagenoise=shocks.wagenoise,
                          shocks.spwagenoise=shocks.spwagenoise,
                          nocomp=nocomp,noeventDI=noeventDI,onlyFT=onlyFT,noflow=noflow,nostockDI = nostockDI, nocontflow = nocontflow,
                          consmoments_simple=moments_simple,
                          workmoments_simple=moments_simple,
                          DImoments_simple=moments_simple)
  
  Jmat<-cbind(Jmat,(moments1-moments0)/(param1-param0))



  write_csv(x=as.data.frame(Jmat),path=paste0("./modeloutput/newguess_jacobian_",modeltype,".csv"),append=FALSE,col_names=FALSE)

  pos<-pos+1
  
}
modeltype<-"simplified"
Jmat<-read.csv(paste0("./modeloutput/newguess_jacobian_",modeltype,".csv"),header=FALSE)

wagemat <- c(1/wageregse^2,1/spwageregse^2,1/wagevarse^2)
consmat <- (1/(consse_simple_all^2))#/sum((1/(consse_all^2)))
workmat <- c((1/(workse_simple^2)),(1/(spworkfullse^2)))
#workmat<-workmat /sum(workmat)
#Including event moments instead of flow and composition moments:
DImat <- c(
  (1/(compDIse_simple^2)),
  (1/(flowDIse_simple^2)),
  (1/(eventDIse_simple^2))
)
weightmat <- diag(c(wagemat,consmat,workmat,DImat))
Jmat<-as.matrix(Jmat)
varcovhat<-as.matrix(varcov)*(1+1/numS)
param_var<-solve(t(Jmat)%*%weightmat%*%as.matrix(Jmat))%*%t(Jmat)%*%weightmat%*%varcovhat%*%weightmat%*%as.matrix(Jmat)%*%solve(t(Jmat)%*%weightmat%*%as.matrix(Jmat))
Lambda<-solve(t(Jmat)%*%weightmat%*%as.matrix(Jmat))%*%t(Jmat)%*%weightmat
rownames(Lambda)<-names(newguess[paramnums])
colnames(Lambda)<-names(moments1)

write_csv(x=as.data.frame(Lambda),path=paste0("./modeloutput/newguess_sensitivity_",modeltype,".csv"),append=FALSE,col_names=FALSE)
write_csv(x=as.data.frame(param_var),path=paste0("./modeloutput/newguess_varcov_",modeltype,".csv"),append=FALSE,col_names=FALSE)
}
